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ABSTRACT 


A  method  is  developed  to  form  beam  spectral  estimates 
from  a  set  of  signals  sampled  from  hydrophone  outputs  of  an 
acoustic  array.  The  signals  must  be  frequency  band  limited 
over  an  interval  [W0,W0+W ]  and  sampled  over  a  time  interval 
T.  A  savings  of  ( WQ/W )  +  1  in  data  storage  and  computation 
time  is  realized  over  conventional  beamforming  analysis 
methods  employing  sampling  rates  equal  to  twice  the  highest 
frequency  in  the  signal  band.  Any  set  of  delays  t,  may  be 
utilized  and  the  beamforming  error  determined  as  long  as 
t,-  <  T/2.  The  computed  beam  spectra  can  be  made  arbitrarily 
close  to  the  actual  beam  spectra  by  choosing  (a)  a  large  num¬ 
ber  of  points  to  be  analyzed,  (b)  small  input  signal  amplitudes 
at  the  end  points  of  the  time  interval  T,  and  (c)  small  actual 
spectrum  amplitudes  near  f=W0  and  f^¥/0+W. 

This  procedure  would  be  particularly  applicable  to  future 
sonar  systems  employing  digital  circuitry  to  beamform  and 
spectrum  analyze  on-line. 


PROBLEM  STATUS 

This  is  an  interim  report  on  the  problem;  work  is 
continuing. 
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DIGITAL  BANDPASS  BEAMFORMING  WITH  THE 
DISCRETE  FOURIER  TRANSFORM 


INTRODUCTION 

The  advent  of  the  Fast  Fourier  Transform  (FFT)  has  greatly  facilitated  computer  calcu¬ 
lations  of  power  spectra.  On-line  spectrum  analysis  of  signals  by  digital  computer  is  expected 
to  play  an  important  role  in  future  sonar  systems.  The  method  presented  here  allows  a  beam- 
formed  spectrum  to  be  computed  from  a  set  of  appropriately  sampled  bandpass  signals  repre¬ 
senting  hydrophone  outputs  of  an  array.  If  the  ratio  of  the  center  frequency  or  carrier  to  the 
signal  band  about  the  carrier  is  much  greater  than  one,  a  conventional  beamforming  analysis 
utilizing  sampling  of  the  data  at  twice  the  highest  frequency  in  the  signal  band  not  only  wastes 
time  and  money  in  providing  many  spectral  estimates  at  frequencies  of  no  interest,  but  requires 
inordinately  large  amounts  of  storage  to  obtain,  a  desired  resolution.  These  high  sampling  rates 
are  normally  used  when  implementing  a  digital  beamforming  analysis  in  order  to  obtain  small, 
discrete  time  intervals  for  delaying  purposes.  A  bandpass  sampling  method  of  beamforming 
would  be  more  efficient,  would  yield  better  resolution  for  the  same  amount  of  input  data 
samples,  and  would  result  in  good  coverage  of  the  useful  bandwidth  utilizing  all  the  available 
spectral  estimates.  Such  a  scheme  is  employed  here  with  the  following  results: 

(a)  The  beamformed  spectrum  is  expressed  as  a  single  equation  utilizing  the  conven¬ 
tional  Discrete  Fourier  Transform  (DFT). 

(b)  The  sampling  rates  are  considerably  lower  than  those  normally  employed  for  digital 
beamforming.  Therefore,  a  substantial  saving  of  data  storage  and  computation  time  results. 

(c)  Sampling  the  signal  at  these  low  rates  does  not  constrain  the  delays  f;  to  be  integral 
multiples  of  the  sampling  interval  and,  in  fact,  any  delay  is  realizable  with  this  method  and  an 
associated  beamforming  error  may  be  calculated  as  long  as  t,-  <  T/2  where  T  is  the  signal 
length  to  be  analyzed. 

(d)  The  use  of  first-order  sampling*  affords  a  very  simple  and  practical  means  of  ob¬ 
taining  the  discrete  data  samples  for  analysis. 


DIGITAL  BANDPASS  SPECTRUM  ANALYSIS 

In  this  section  we  review  the  concept  of  bandpass  sampling  and  develop  a  method  for 
estimating  the  power  spectrum  of  a  bandpass  function  g(t).  The  results  of  this  section  are 
then  used  in  the  next  section  to  develop  and  construct  the  beam-spectrum  equation. 


*A  sampling  is  defined  as  first-order  if  the  sample  points  are  equispaced.  See,  for  example,  Ref.  1. 
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Consider  a  function g(t)  which  is  band  limited  to  the  frequency  interval  [W0,W0+W] . 
The  function  g(t)  may  be  reconstructed  exactly  (2)  from  an  infinite  set  of  first-order  samples 
taken  at  a  rate  2 W  only  if  WQ  =  CW,  where  C  =  1,2.3,... .  Although  higher  order  sampling 
permits  reconstruction  of  the  original  signal  for  less  severe  constraints  on  W0  and  W,  we  will 
only  consider  the  case  for  which 


W0  =  CW,C  =  1,2,3,,.. 


(1) 


since  first-order  sampling  is  more  amenable  to  current  operational  digitizing  systems.  The 
complexity  of  higher  order  sampling  schemes  makes  first-order  sampling  a  very  simple  and 
practical  means  of  obtaining  the  discrete  data  samples.  In  a  practical  situation  for  which  WQ 
is  not  an  integer  multiple  of  the  signal  bandwidth,  W  can  be  easily  chosen  slightly  larger  than 
the  actual  signal  bandwidth  so  that  Eq.  (1)  is  satisfied. 


Let  samples  of  g{t)  be  taken  at  a  rate  2 W  and  let  WQ  =  CW,  where  C  is  a  positive  integer. 
The  conventional  Nyquist  rate  for  sampling  a  function  which  is  band  limited  to  the  interval 
[0,Wo+W]  is  equal  to  2{W0+W).  Thus,  by  sampling  at  a  rate  2 W  (“submultiple”  or  “band¬ 
pass”  sampling),  a  factor  of  (WQ/W)+1  =  C+l  savings  in  data  storage  is  realized.  We  define 
the  reconstruction  gR  (t)  of  g(t)  as  in  (2)  by 


where 


2(W0+W)  sine  2 (Wc+W) 


2 WQ  sine 


sincX  = 


sin  nX 
7 tX  ' 


(2) 


It  can  be  shown  that^(t)  is  equal  to  g(t),  exactly,  for  all  t. 

In  addition  to  being  band  limited,  let  g(t)  be  time  limited  to  the  interval  [0,T] .  There 
are  now  2 WT  =  N  samples  obtained  from  sampling  ^(f)  in  the  interval  [0,T] .  The  reconstruc¬ 
tion,  Eq.  2,  for  #(f)  remains  the  same,  except  that  the  limits  on  the  sum  are  changed  to 
n  =  0,  and  n  =  (N—l).  Even  though  a  function  cannot  strictly  be  both  band  limited  and  time 
limited,  these  conditions  may  be  met  approximately  (3).  From  this  point  we  will  consider 
only  those  g(t)  whose  values  are  negligible  outside  some  time  interval  of  length  T  and  whose 
spectrum  values  are  also  negligible  outside  the  frequency  interval  [W0,W0+W } .  The  function 
gR  (t )  is  now  approximately  equal  to  g(t)  over  the  interval  T.  The  closeness  of  the  approxi¬ 
mation  depends  on  the  smallness  of  the  spectrum  of  g(t)  outside  the  interval  [W0,W0+W] 
and  the  smallness  of  g(t)  outside  T  (see  App.  A). 
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Define  the  Fourier  transform  pair  to  be 


oo 

ng(t))=G(f)=  J  g(t]e-W*dt 

—oo 

oo 

ri[G(/0]=^(f)=  |  G(f)eW*df 


where  t p  denotes  the  transform  operator  and  f  is  measured  in  Hertz.  Two  transforms  which 
will  be  used  in  the  development  are 


i^[sinc(Zf)]  = 


— ,  \f\  <— 
X  1  2 


0,  \f\>j 


^[^(f+or)]  =c2  nif<*$[g(t)]  =  e™f«G(f). 

We  now  proceed  to  derive  an  expression  for  the  transform  of  the  bandpass  sampled  .  -. 
time  function  g(i)  in  terms  of  the  Discrete  Fourier  Transform  (DFT).  First,  we  apply  the 
operator  $  to  gR  (f )  to  obtain 


GR(f)  =  ^[gR(t)]  = 


N-l  .  >1 

SE  *  (is)  <r2,r'w2w  )  ,  Wo  <  \f\ < 

n=0  '  '  •  >  .  ( 


elsewhere 


Equation  (4)  expresses  the  transform  of  gR  (f)  as  a  function  of  continuous  frequency  f. 
However,  the  computer  operates  on  a  discrete  set  of  data,  i.e.,  on  the  values  g(n/2W)  for 
n  =  The  DFT  resulting  from  the  operation  consists  of  N  output  values  at  the 

frequencies 

f  =*.=  2 W 
fk  T  N  k 


where  k  is  an  integer  satisfying  the  relation  W0T  <  jfe|  <  (W0+IP;T. 
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Substitution  of  fk  for  fin  Eq.  (4)  yields 


c*(4)- 


W0<\fk\<W0+W  I 


^  0,  elsewhere 


(5) 


where  fk  =  k/T. 

The  sufficiency  conditions  that  enable  the  function  GR  (4 )  to  be  made  arbitrarily  close 
to  the  “correct”  transform  G(f),  evaluated  a.tf~fk,axe  shown  in  App.  A.  These  conditions 
show  the  amount  of  data,  the  time  domain  weighting,  and  the  degree  of  bandpass  filtering 
that  is  necessary  to  meet  a  desired  degree  of  accuracy  in  GR  (fk ). 

The  form  of  Eq.  (5)  looks  much  like  the  DFT.  However,  the  DFT  is  usually  computed 
for  ft  =  0,1.2,. The  bounds  on  the  index  ft  in  Eq.  (5)  are 


W0T<\k\<(W0+W)T. 
Since  2WT  =  N  and  WQ  =  CW,  we  have 


(c*m 

O  '  '  O 


(6) 


Therefore,  adjustments  must  be  made  in  order  to  properly  evaluate  the  spectrum,  of  g(t)  by 
means  of  the  conventional  DFT. 


Since  C  is  an  integer,  we  can  write 


C  =  2m  —  £ 


where  m  is  a  positive  integer  and  £  =  0orlasCis  even  or  odd. 


If  we  introduce  a  new  index  ft'  such  that 


'  (CN/2)  +  ft'  =  mN  —  (RN/2)  +  ft',  ft'  =  0,1,...,  {N/2)-l 
I  —{CN/2)  +  ft'  -  -mN  +  {m/2)  +  ft',  ft'  =  -A72,...,-2,-l 


then 


for  0<|fc’|«W2 


Substituting  Eq.  (7)  into  Eq.  (5)  and  noting  that  e  27!ik  n-,N  is  periodic  in  ft'  with  period  N, 
and  that  e~nin *  -  e*nini,  we  have 
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(jf  £  (g  fjf)  e^e-2 nink’/N)  ,  k'  =  0,1,-, N-l 

vr(/*)  =  -S  n-0  \  ' 


all  other  k' 


where 


CN  V 
2T  +  T’ 


k'  =  0,1,2 . (Nf  2)-l 


_M+|,  v-w/a . jh 

2  T  T  ’ 


In  order  to  reduce  Eq.  (8)  to  the  form  of  the  DFT  we  must  consider  the  two  cases  for  C 
even  or  odd.  When  2  =  0  (C  is  even)  Eq.  (8)  reduces  to 

j?Z  (g(j?)  ^n,nk' ’/N)  >  »'  -  04 . N-l 

gr(4)  =  "s  -o  '  V  /  ^ 


all  other  fe' 


where 


miV+fe' 


,  &'  -  0,1,...,(AV2)-1 


=£££.  v-w/wM 


This  is  simply  the  equation  for  the  DFT  of  g(nT/lV),  so  that 

[«(^)1  ,4'  =  0,1,...,JV-1.  (10) 

Thus  for  C  an  even  integer,  the  “correct”  Fourier  transform  of  the  bandpass  function 
5(f)  sampled  over  the  interval  T  may  be  approximated  by  means  of  the  DFT  of  the  N  band¬ 
pass  samples  of  5(f).  Figure  1  shows  graphically  a  spectrum  in  the  three  stages  of  processing: 
(a)  the  original  spectrum  G(f)  for  C  an  even  integer;  (b)  the  spectrum  of  G(f)  after  sampling 
at  a  rate  2 W;  and  (c)  the  spectrum  obtained  from  the  DFT  of  the  N  sampled  values  of  g(t). 


RELATIVE  AMPLITUDE 
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When  £  =  1  (C  is  odd)  we  find* 


(o) 


(b) 


Fig.  1— The  ideal  spectrum  G{f)  for  an  even  value  of 
C  is  illustrated  in  (a).  After  sampling  at  a  rate  2W, 
the  spectrum  becomes  as  shown  in  (b).  Finally,  the 
spectrum  obtained  from  the  Discrete  Fourier  Trans¬ 
form  (DFT)  of  the  N  sampled  values  oig(t)  is  shown 
in  (c).  These  three  sketches  show  the  amplitude  spec¬ 
trum  in  three  stages  of  processing. 


Thus  for  C  odd,  the  data  may  be  spectrum  analyzed  using  the  DFT  if  the  sign  of  every 
other  data  point  is  first  changed.  If  the  percentage  bandwidth  of  the  signal  is  small  enough, 
it  may  be  advantageous  to  increase  W  slightly  to  make  C  even.  Then  the  (~l)n  term  can  be 
avoided. 


♦The  factor  enin  =  (-l)”  is  equal  to  where  f  =  nT/N  and  W  =  N/2T,  corresponding  in  the  frequency 

domain  to  a  shift  of  W;  i.e., 


ip[e2nitWg(t)]  =»  G(f~W). 
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The  Inverse  Discrete  Fourier  Transform  (IDFT)  of  GR  Ofe )  for  C  even  yields  the  original 
time  series  g(nT/N).  The  function  gR  (t)  may  then  be  formed  by  using  the  appropriate 
sampling  interpolation  formula.  The  IDFT  of  GR(f ^ )  for  C  odd  yields  the  time  series 
(-l)n  g(nT/N).  The  function  gR  (t)  may  be  formed  by  first  changing  the  sign  of  every  other 
data  value  and  then  applying  the  sampling  interpolation  formula. 

The  power  density  spectrum  S(fk )  is  formed  from  the  Fourier  transform  of  g(t)  by 
letting 


Q/X  \  —  _  r*  Sr  \ 

•-(4)  -  T  ~RUkt^n  Wfe  f 


where  f  denotes  the  complex  conjugate. 

Thus,  in  this  section  we  found  that  the  Fourier  transform  of  a  function  g(t),  appropri¬ 
ately  time-  and  frequency-band  limited,  can  be  computed  from  a  finite  set  of  bandpass 
sampled  data  of  g(t)  by  means  of  the  DFT,  and  at  a  computer  storage  savings  proportional 
to  the  ratio  of  the  highest  frequency  in  the  signal  band  to  the  bandwidth  of  the  signal.  Ji 
the  number  of  data  samples  N  is  a  power  of  two,  the  FFT  may  be  employed  to  calculate  the 
DFT. 

It  now  remains  to  incorporate  an  arbitrary  delay  into  g(t)  for  the  purpose  of  beam¬ 
forming.  This  will  be  done  in  the  next  section. 


DIGITAL  BANDPASS  BEAMFORMING 

The  previous  section  presented  and  developed  some  properties  of  bandpass  functions 
relating  to  the  DFT.  Once  these  methods  are  at  hand,  it  is  a  fairly  simple  matter  to  incorpo¬ 
rate  a  time  delay  tQ  into  the  equations  and  develop  an  expression  for  the  transform  of  the 
delayed  function  g(t-tD).  After  doing  this,  the  expression  will  be  generalized  to  include  a 
set  of  signals,  each  with  its  own  delay,  representing  the  outputs  of  hydrophones  in  an  array 
to  be  beam  formed  and  spectrum  analyzed. 

Basically,  a  beam  is  formed  from  signals  which  are  delayed  by  predetermined  values  and 
then  summed.  These  delays  can  be  introduced  in  either  the  time  domain  or  the  transform 
(frequency)  domain  of  the  signal  because  of  the  relation 


where 


(12) 


G(f)=+[g(t)]. 
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Introduction  of  a  time  delay  tQ  in  Eq.  (4)  results  in 


CdR(f)  =  'P[gR(t-t0)]  = 


^  _  e~2irift0 
2W 


£  fg  (jf)  e~2^fn/2w\  ,  WQ  <  |/1  <  W0+W^{ 


0,  elsewhere 


where  G^  (f)  is  the  transform  of  the  delayed  reconstruction  function. 
Now  let  f  -  fu  =  k/T.  as  in  the  Drevious  section.  Then. 


f  Z  e~2m,0k/T  ^  (gl*f\  «-»»/•»] 

«-0  V  '  '  ' 


(C+1)N 

2 


v  0,  all  other  k  J 

It  is  clear  from  Eq.  (14)  that  we  are  forming  the  transform  GR  (fk )  of  the  delayed  recon¬ 
struction  function  from  sampled  values  of  g(t)  for  0  <  t  <  T.  However,  the  nature  of  the 
sampling  process  is  such  that  g(nT/N)  appears  to  be  cyclic  in  n  with  period  N,  and  GR(fk’)is 
cyclic  in  k'  (see  previous  section)  with  period  JV.  A  delay  of  tQ  will  therefore  cause  the 
function  g(t),  0  <  t  <  T,  to  be  delayed  into  the  new  function  (see  Fig.  2) 


**(0- 


^(0  = 


ta<t<T 


^g(T~t0+t),  0<t<to 
rg(t+|fJ),  0<t<(T-\to\) 


g(t—T+\t0\),  (T-\t0\)<t<T 


for  tQ  >  0 


for  tQ  <  0 


Therefore,  those  values  of  the  delayed  reconstruction  function  gR  (tn-tQ )  computed  for  the 
interval  [0,fo]  for  tQ  >  0,  or  for  the  interval  [7Ht0  l,T)  for  t0  <  0,  may  not  be  accurate 
since  this  part  of  the  function  has  been  “folded  over”  or  “cycled”  from  the  other  end  of 
g(t).  However,  if  the  quantity  tD  is  much  less  than  T,  any  error  in  the  delayed  function 
gR(tn-t0)  on  the  interval  [0,|fo|]  or  [T-]t0|,T]  will  detract  negligibly  from  the  correctness 
of  the  function  on  the  interval  [0,T]  (see  App.  B).  In  fact,  for  most  practical  acoustic 
underwater  array  problems,  tp  will  be  on  the  order  of  milliseconds  and  T  will  be  several 
seconds,  so  that  (f0/T)  1CT  2.  An  error  equation  is  developed  in  App.  B  which  includes 
the  effects  from  “cyclic”  delaying.  It  is  shown  in  App.  C  that  in  at  least  one  specific  case, 
the  cyclic  delaying  error  is  small  for  delays  close  to  T/2. 


To  reduce  Eq.  (14)  to  the  form  of  the  DFT,  we  proceed  as  in  the  previous  section. 
Since  C  is  an  integer  we  can  write 


NRL  REPORT  7359 


9 


Fig.  2— A  delay  of  tQ  will  cause  the  time-limited  function  g(t) 
[sketch  (a)]  to  be  delayed  into  the  new  function  for 

tQ  >  0  [sketch  (b) ]  or  for  tQ  <  0  [sketch  (c)]. 


C  =  2m— £ 

where  m  is  a  positive  integer  and  fi  =  0  or  1  as  C  is  even  or  odd.  Introduce  a  new  index  ft' 
such  that 

qn 

mN  +  ft',  ft'  =  0,l,2,...(JV/2)-l 
z 

CAT 

-i mN  +  ft',  ft'  =  -N/ 2,. ..,-2,-1 

z 

We  substitute  ft'  for  ft  in  Eq.  (14)  and  then  specialize  to  the  two  cases  of  even  or  odd  C. 
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For  £  =  0  (C  even)  we  find 

leJTrJ>(mN*k  ’  DPT*'  [  .  »'  =  0.1 . (W2M 

aRVk)  =  *  '  DFTfc.  f  ,  k'  =  N/2,...Jt-l  ►  ' 

0,  all  oilier  ft' 

^  ^  ttfil 

%  *  / 

Thus  for  C  even,  the  transform  of  the  delayed  function g(t~t0)  is  approximated  by  com¬ 
puting  the  DFT  of  the  set  of  values  g(nT/N)  and  then  multiplying  the  fe'th  value  by  the  ap¬ 
propriate  exponential  phase  factor. 

For  2  =  1  (C  odd)  we  have 

■f  «  T  °  <m‘V  "  <Ar/2)  *  ‘  ’  DFT..  [(-!)»  ,  ft'  =  0 . (MH 

k..m . ^  - 

0,  all  other  ft' 


Thus  for  C  odd,  the  transform  of  the  delayed  function  g{t~tQ)  is  approximated  by  computing 
the  DFT  of  the  set  of  values  (-1)*  g(nT/N)  and  multiplying  the  fe'th  value  by  the  appropriate 
exponential  phase  factor. 


Equation  (16)  or  (17)  allows  us  to  approximate  the  Fourier  transform  of  a  delayed 
time  function g(t-tQ)  by  operating  on  bandpass  samples  of  g(t)  taken  from  t  =  0  to  t  ~  T. 
The  only  new  restriction  is  that  the  delay  tQ  should  be  small  compared  to  the  length  of 
signal  T  to  be  analyzed. 


In  order  to  generalize  the  equations  for  the  beamforming  case,  consider  an  array  of  p 
hydrophones,  in  which  the  rth  hydrophone  receives  a  signal  f{nT/N ),  r  =  1,2,...  jj.  Let  the 
delay  to  be  introduced  into  the  rth  received  signal  be  ir  and  let  (tr)max/T  be  much  less  than 
unity.  For  C  an  even  integer  the  transform  of  the  formed  beam  (computed  from  all  receiver 
output  signals  1  to  p,  each  having  been  sampled  at  a  rate  2  W  from  0  to  T)  is  approximated 
by 
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C  T  P  l2wifr  , 

'e  T  (mN*k)BVTk. 


jfL 


r-1 


Gl0am(4H  T  p  \  -™<r  , 

lv  L 

r-l 


e 


Y4?\ 

\iV  / 


k’  =  0,1 . (iV/2) — 1 


,  k’  =N/2t...,N-l 


> 


0,  all  other  k1 


ft  Q\ 
V^w; 


This  is  the  final  result  of  the  analysis.  The  beam  spectrum  is  readily  calculated  by  means 
of  Eq.  18  which  utilizes  the  DFT  and  the  bandpass  sampled  date  from  the  individual  hydro¬ 
phones  in  the  array.  Equation  (18)  couid  be  implemented  with  a  digital  computer  employ¬ 
ing  the  Fast  Fourier  Transform  to  compute  the  DFT  of  the  input  signals.  To  save  even  more 
computation  time,  only  those  Fourier  estimates  at  frequencies  fk '  of  interest  need  be  multi¬ 
plied  by  the  exponential  phase  factors  and  summed.  A  wide  choice  of  analysis  ‘‘windows”  is 
available  (3),  and  incorporation  of  such  a  function  in  Eq.  (18)  would  serve  the  dual  purpose 
of  reducing  spectral  leakage  and  providing  the  necessary  time-limiting  approximations  on 
g(t). 


SUMMARY 

Signals  which  are  band  limited  to  the  frequency  interval  [W0,W0+W\  can  be  sampled 
from  f  =  0tof  =  Tata  rate  2  W  and  digitally  processed  according  to  Eq.  (18)  to  obtain  a 
beamformed  spectrum  of  the  inputs.  A  savings  of  (W0/W)+l  in  date  storage  or  an  improved 
frequency  resolution  is  realized  over  the  conventional  digital  beamforming  methods  because 
of  bandpass  sampling.  The  computed  spectrum  can  be  made  arbitrarily  close  to  the  actual 
spectrum  by  choosing  a  large  number  of  points  analyzed,  small  input  signal  amplitudes  near 
t  -  0  and  t  =  T,  and  small  spectrum  values  near  f=WQ  and  f  -  Wa+W.  The  delays  are  not 
constrained  to  be  integral  multiples  of  the  sampling  interval  At.  Any  set  of  delays  t,-  may 
be  used  and  a  beamforming  error  term  calculated  as  long  as  t,-  <  T/2.  This  method  could 
be  used  for  digital  on-line  analysis  of  narrow-band  signals  constituting  outputs  of  a  receiv¬ 
ing  array  in  a  sonar  system,  or  it  might  be  used  as  a  time-  and  space-saving  date  processing 
technique  employing  a  conventional  digital  computer. 
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APPENDIX  A 


BANDPASS  SPECTRUM  ERROR 


There  are  five  sources  of  error  which  contribute  to  the  total  error  in  the  beamformed 
amplitude  spectrum  ).  Only  one  of  these  errors  originates  from  the  delay  properties 

of  the  processing.  The  rest  are  present  in  the  Discrete  Fourier  Transform  (DFT)  regardless 
of  any  delay  introduced.  The  five  errors  are 


t  i .  •  _ j_  l  _  _______  m  /<f  /jp\  •  __  _  x  ,  f, '  i  t  t  .  it*  *  *  *  » 

spcCuulu  aliasing  cuui,  picacni  ucOaUoc  Oi  Uy  }  is  iiOi>  S&iCUy  L'iUIUliUllUJU  lO  nr , 


b.  truncation  error  from  the  time  domain  caused  by  sampling  g(t)  only  over  the 
interval  T, 

c.  spectral  leakage  through  sidelobes  of  the  analyzing  window,* 

d.  error  from  spectrum  curvature,  and 

e.  a  time  domain  cycling  error  caused  by  the  cyclic  delaying  properties  of  the  beam¬ 
forming  process. 


In  this  appendix  we  discuss  the  first  four  errors,  which  contribute  to  the  total  error  in 
the  bandpass  spectrum  computed  via  the  DFT.  (In  App.  B  we  will  discuss  the  error  in  the 
beamformed  bandpass  spectrum.)  The  errors  from  leakage  and  curvature  are  discussed 
qualitatively;  the  aliasing  and  truncation  errors  will  then  be  accounted  for  quantitatively  by 
combining  their  effects  into  a  single  error  formula. 

The  degree  of  spectral  leakage  that  can  be  tolerated  is  largely  dependent  on  the  char¬ 
acteristic  of  the  spectra  to  be  analyzed.  If  the  spectra  fire  relatively  flat  in  the  band  of  inter¬ 
est,  leakage  of  energy  through  sidelobes  of  the  analyzing  spectral  window  would  not  be  of 
concern.  But  for  spectra  with  a  large  dynamic  range  of  amplitudes,  or  for  one  with  sharp 
separated  spikes,  leakage  can  result  in  false  estimates  of  power.  If  leakage  is  a  concern,  the 
Hanning  window  is  an  excellent  smoothing  window  which  may  be  applied  either  as  a  multi¬ 
plication  of  the  time  data  or  as  a  convolution  of  the  complex  Fourier  amplitude  spectrum  (3). 
Although  the  resolution  is  reduced,  the  peak  level  of  the  Hanning  spectral  window  at  the 
first  sidelobe  is  more  than  30  dB  down  from  the  peak  level  at  the  main  lobe. 

A  restriction  or  the  curvature  of  the  spectrum  is  made  by  requiring  that  the  energy  in 
the  ideal  spectrum  be  almost  constant  over  the  analyzing  window  width.  Here  “analyzing 


♦The  analyzing  window  may  be  thought  of  as  a  spectral  weighting  function  of  unity  gain  which  is  convolved 
with  the  ideal  spectrum  G(f)  to  obtain  the  estimates  of  spectral  amplitude. 
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window  width”  refers  to  the  width  of  the  main  lobe  of  the  spectral  analyzing  window.  This 
requirement  assures  that  an  unbiased  estimate  of  the  true  spectrum  will  be  obtained  (Ref.  3, 
p.  8).  A  user  must  choose  an  analysis  time  T  long  enough  to  assure  enough  resolution  or 
minimal  window  width. 


Now,  to  obtain  the  quantitative  estimate  of  the  error  from  aliasing  and  truncation,  we 
follow  the  procedure  developed  by  Moseley.  Moseley  calculated  these  errors  in  the  DFT  of  a 
time  function  g(t )  sampled  over  the  interval  T.  The  spectrum  of  g(t)  in  his  development  is 
low  pass  and  is  bandlimited  to  a  frequency  W  (Ref.  3,  App.  B).  We  will  use  the  results  of 
his  analysis  to  obtain  a  measure  of  the  aliasing  error  and  then  to  obtain  a  measure  of  the 
truncation  error  for  our  bandpass  spectrum.  Tn  order  to  apply  his  results  to  the  bandpass 
case  we  proceed  to  split  the  spectrum  into  two  parts.  This  enables  us  to  compute  the  error 
in  the  DFT  caused  by  aliasing.  Let 


n/*\  —  a'{f\ 
^  (i ;  Ki ; 


>  G"(f) 


(Ai) 


where 


and 


Recall  that 


G'(f)  = 


w0<f<w0+w' 

elsewhere 


G"(f) 


f  G(f), 


-w0-w<f<-woj 

elsewhere  j 


G(f)  = 


G(f),  W0<\f\<(W0+W) 
<  0S  elsewhere 


We  proceed  to  find  the  aliasing  error  in  G'{f)  and  G"(f)  and  then  add  these  together  to  ob¬ 
tain  the  error  in  G(f). 


Consider  a  frequency  translation  of  G'(f)  to  G'(f),  where  f'  =  f~W0-  (W/2)*.  Then 


G'(f')  = 


'  G\f  \  0<\f'\<W/2  ' 

"  ?•. 
«  0,  W/2  <  |/"'|  <oo 


(A2) 


*This  frequency  shift  merely  amounts  to  a  complex  phase  m  .'tipli cation  of  the  time  data. 
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Assume  that  | G'(f')\  approaches  zero  at  least  as  rapidly  as 

Jrz  for  (A3) 

1 1  I 

where  a  >  1.0  and  e  is  some  small  constant.  Denoting  the  aLiasing  error  by  7i ,  Moseley 
calculates  this  error  in  G'(f')  to  be  given  by 


iTil  < 


\e\2«+H{ct) 

(SF)a 


(A4) 


where  SF  is  the  <wmr»lintr  frenueney  (2!^  for  OUT  CSSS),  Slid  j(s)  IS  the  RlCuiuTifi  ZctS  fUIiCtlGn 
defined  by 


oo 

w_.\  _  \  '  n— n> 

-  Z_*  x  • 


Similarly,  after  translating  G"(f)  to  G"(f"),  where  f"  =  f+W0+(W/2),  we  find  that  the  aliasing 
error  in  G"(f")  is  also  less  then  y1 .  This  time  we  require  that  |G"(f  ")|  approach  zero  at  least 
as  rapidly  as 


tor  W/2  <  in  <  °°- 

For  our  bandpass  spectrum,  the  restrictions  on  G'{f')  and  G"(f")  require  that  G{f)  approach 
zero  at  least  as  rapidly  as 


llfl  -  W0-(W/2)|  “  for  lfl  outside  lw°’w°+w]  ■ 

Then  these  errors  may  be  added  to  obtain  an  approximate  estimate  of  the  aliasing  error  in 
the  bandpass  spectrum: 


l7ipl  < 


|e!2a+2r(a) 

(SF)a 


(A5) 


This  error  is  twice  the  error  in  the  low-pass  spectrum  bandlimited  to  W  because  for  the  band¬ 
pass  case  we  have  two  frequency  cutoff  regions,  while  for  the  low-pass  case  there  is  only 
one.  While  there  are  only  C  aliased  bands  from  frequencies  less  than  W0,  and  an  infinite 
number  of  aliased  bands  from  frequencies  above  Wq+W,  the  falloff  of  the  out-of-band  energy 
is  usually  fast  enough  so  that  doubling  the  low-pass  aliasing  error  is  a  good  approximation 
to  the  aliasing  bandpass  error.  In  any  case  the  estimated  bandpass  aliasing  error  is  a  valid 
upper  bound  of  the  actual  aliasing  error. 

To  obtain  an  estimate  of  the  truncation  error,  we  require  thatg(f)  approach  zero  at 
least  as  rapidly  as  vi\t\$  for  t  outside  the  sampled  interval  T,  where  |3>1.0  and  v  is  some 
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small  constant.  Mosely  calculates  the  truncation  error,  denoted  here  by  72,  to  be 
given  by 


l72l< 


\p\2? 

(/KLKT-A^-1 


(A6) 


where  A  =  1/SF  is  the  interval  between  samples. 


The  total  aliasing  and  truncation  error  in  the  DFT  of  the  bandpass  function  g(t)  is  then, 


if  m  cmi  -r  |f|2°*2»a> , 

Gr(-  l-  {SF)<*  (0-l)(T-A)P_1 ' 


The  aliasing  error  may  be  reduced  by  choosing  a  larger  sampling  frequency  SF,  although  in 
the  bandpass  case  it  must  be  an  integral  multiple  of  2 W.  The  aliasing  error  could  be  further 
reduced  if  the  data  were  filtered  with  a  filter  of  steeper  slope  such  that  the  signal  in  the  band¬ 
width  of  interest  was  not  affected.  In  this  case  the  quantity  ft  would  be  larger.  The  trunca¬ 
tion  error  may  be  reduced  by  choosing  a  longer  analysis  time  T. 


APPENDIX  B 


BEAMFORMED  SPECTRUM  ERROR 


An  estimate  of  the  aliasing  and  truncation  error  present  in  the  DFT  of  a  bandpass  func- 
g(t)  was  determined  in  App.  A.  here  we  examine  the  error  in  the  beamformed  spectrum 

Consider  a  set  of  band-limited  functions  ^(t)  for  r  =  l,2,...,p,  all  of  which  are  identical 
after  time  shifting  by  fr,  where  we  let  t1  =  0  be  the  reference  channel  chosen  so  that  the 
remaining  ir  are  positive.  Then, 


L 


0? It — t  ^  —  nrrl  I 

5  ~  /  X~ Cf  /• 


If  the  set  Gr(f)  are  the  Fourier  transforms  of  the  we  have 


P 

£  (e-2 =pGltf). 
r«»l 


(Bl) 


(B2) 


We  now  sample  each  gr( t )  in  the  interval  [— T/2,T/2]  at  a  rate  2 W  =  N/T  =  1/A  and  beam- 
form  via  the  DFT  using  the  method  described  in  the  main  body  of  this  report.  Let  GrR  if)  be 
the  set  of  spectra  computed  by  the  DFT’s  of  the  data  before  the  beamforming  is  carried  out. 
Then, 


P 

Gflam(n=  ]T  [e-WtrG^f)]  and 

r<=  1 


<r>-L 


[e-2KiftrGr(f)] 


are  the  beam  spectrum  computed  via  the  DFT  and  the  ideal  beam  spectrum,  respectively. 
Using  Eq.  (B2)  and  defining 

P 

£  spGjj!0.  ns) 

r»  1 

the  difference  between  G^eam  (f)  and  G1*5®"1  (f)  is  pG^  (f)  -pG1  (f).  The  factors  of  p  are  in¬ 
dicative  of  the  signal  gain  that  we  expect  from  the  beamforming  process. 

Moseley  derives  an  equation  for  the  truncation  and  aliasing  error  present  in  the  DFT 
of  a  sampled  low-pass  time  function  gn{i).  His  error  equation  is  expressed  as  a  difference 
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between  the  spectrum  GR  (f)  computed  via  the  DFT  and  the  ideal  continuous  spectrum 
G(f),  where 


G(f)  = 


oo 

— oo 


and 


(N/2)~l 

GR(f)  =  teiAnJe-atri/A*] 

n=~N/2 

with  A  =  1/SF  ~  T/N,  and  SF  is  the  sampling  frequency. 

A  brief  explanation  of  his  method  follows.  By  sampling  g(t)  with  an  infinite  Dirac 
comb  function,  and  taking  the  continuous  transform,  he  obtains  for  \f\  <  1/(2 A),  the  true 
transform  G(f)  plus  an  infinite  number  of  aliased  spectral  terms.  Then  he  takes  into  account 
time  truncation  by  subtracting  samples  of  g(An)  for  n  =  N/2  to  °°  and  obtains  a  series  of 
error  terms  caused  by  truncation.  Thus  he  obtains  a  measure  of  the  aliasing  and  truncation 
error  between  GR  (f)  and  G(f). 

When  we  apply  our  beamforming  procedure  to  this  difference,  we  obtain 


£  j  e~2niftr[GR (f)  -  Gr(f)] }  --Gg™(n-Ghe*m(f)=pG1R(f)-pGl(f). 


r=  1 


Thus  the  exponential  multiplication  and  summing  of  the  Moseley  difference  errors  will  give 
us  a  measure  of  the  error  between  the  DFT  and  ideal  beam  spectra.  Doing  this,  we  obtain 


GbCam(n_c?boam(/)=  £  J  £  e~2niftr  lf-jj  +  Gr 

r=l  1,2=1  L 


-A  gr(nA)e-2«if(nA+tr) 

r=l  I 

l,  n=N/2 


~{N/ 2)-l 


gr(nA]e~27tiftn&+tr) 


> 


(B4) 
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The  right-hand  side  of  Eq.  (B4)  consists  of  two  major  terms  which  we’ll  denote  as  7j  and 
72  separately  consider  by  virtue  of 

I <?£“"(/)-  G^(f)\  =  l7i  +  72i  <  l7i !  +  l72i-  (B5) 

The  first  term,  further  reduced  in  Eq.  (B6)  below,  represents  the  contribution  to  the  beam¬ 
forming  error  caused  by  spectrum  aliasing: 

Note  that,  because  the  multiplicative  phase  terms  are  simply  phase  shifts  not  affecting  the 
amplitudes  of  the  spectra,  our  upper  bound  estimate  of  the  aliasing  error  is  not  affected  at 
all  by  the  introduction  of  the  delays  tr. 

As  in  App.  A,  assume  that  the  Gr(f)  each  approach  zero  at  least  as  rapidly  as 


— -----  for  If  I  outsidp 

\\f\-W0-{W/2)\  - -  0  J 

Then,  except  for  a  factor  of  p.  the  quantity  7j  is  simply  the  aliasing  error  previously  defined 
in  App.  A  and  is 


l7i!< 


(SF)“ 


(B7) 


Recall  that  SF  is  equal  to  2W  for  the  bandpass  case,  and  f (a)  is  the  Riemann  zeta  function. 


The  second  term  72  includes  the  effects  of  time  domain  truncation  and  cycling  error. 
Let  n1  =  n  +  trfA.  Then  we  have 
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(BIO) 
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The  first  term  on  the  right-hand  side  of  Eq.  (BIO)  is  p  times  Moseley’s  truncation  error.  It 
is  less  than  the  quantity 


p\v\  & 

(/j-ixT-Ay*-1 


(Bll) 


for  g(t)  <  where  t  is  outside  the  interval  [-772,772] ,  (3  >  1.0,  and  v  is  some  small 
constant  (see  App.  A). 

The  second  term,  hereafter  denoted  by  7^,  is  caused  by  the  cyclic  representation  of 
gn(t).  When  a  delay  of  g{nA)  is  effected,  a  “cycling”  takes  place,  resulting  in  a  time  inter¬ 
val  where  the  delayed  function  is  inexact.  For  example,  for  positive  tr,  information  in  the 
interval  [T/2-fr,T/2]  is  “lifted”  and  “inserted”  into  the  interval  [-T/2,-T/2+tr] .  For  small 
values  of  delay  tr,  this  error  interval  is  small,  as  is  the  error  itself.  This  can  be  seen  in  Eq. 
(B8)  by  noting  that  if  tr  =  0,  no  additional  terms  are  present  to  contribute  to  y\\  if  tt  =  A, 
only  one  additional  term  contributes  from  each  sum,  etc. 

In  order  to  reduce  and  evaluate  7|,  note  that 


Now  letgr(t— tr)  approach  zero  at  least  as  rapidly  as 

for  T/2  <  |f|  <  00 

where  0  >  1.0,  and  v  is  some  small  constant.  Recalling  that  gr(t)  approaches  zero  for 
t  <  —T/2,  this  is  equivalent  to  letting  gr(t)  approach  zero  for  t  outside  the  interval 
[~T/2,T/2— tr]  for  positive  fr.  In  order  to  maximize  the  error  estimate,  choose  ;max 
equal  to  the  largest  delay  from  the  set  fr.  Then, 
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Finally,  since 


(AT/2)-l+fmax/A 


(|nl|)-0<^« 


nl=!V/2 


(B14) 


/i_i  iy-g 

vi«- 1;  y 


nl=(A72)+l-tmax/A 


A 


(N  t~~\-P 


we  arrive  at  the  expression 


■  L+  _ 1 _ 

(JL'f  1  (1+2_WV 

\  2  /  \  \  N  (T/2)j 


^15101 


where  we  have  used  the  relation  NA  -  T.  It  can  be  seen  from  Eq.  (B15)  that  the  cycling  error 
may  be  decreased  by  choosing  the  analysis  time  T  to  be  large. 

The  estimate  of  cyclic  delaying  error  developed  here  is  good  for  all  practical  situations 
where  the  maximum  delay  tm  ax  is  much  less  than  the  total  length  T  of  signal  being  analyzed.* 
In  fact  the  model  is  good  for  delays  up  to  half  the  analysis  interval  (see  App.  C),  an  extreme 
delay  situation  rarely,  if  ever,  encountered  with  real  arrays  and  signals.  At  tmax  =  (T/2)  +  A, 
the  error  term  I72I  becomes  infinite.  This  is  because,  for  this  particular  delay,  we  require 
that  g(t)  be  negligible  outside  the  interval  [-T/2,0] .  The  function  v/\t\&  has  a  singularity  at 
t  =  0,  however,  thus  preventing  negligibility  of  g{t). 

We  may  now  add  together  the  three  error  terms  (aliasing,  truncation,  and  cycling)  to 
obtain 


IG&eam(fl  "  I  <Pkl(2gp)lf(a) 


(p-I)(p-A)M-i  (77 2f 


Vi  + 

1 

\ 


1  \ 
i  2  i  yf 


(Bis) 


*  Recall  from  the  section  on  Digital  Bandpass  Beamforming  that  a  ratio  of  tmax/T  on  the  order  of  10  2 
was  considered  typical. 
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Thus,  the  error  in  the  beamformed  spectrum  is  the  sum  of  three  error  terms:  the  first  two 
are  encountered  in  the  bandpass  spectrum  analysis,  with  or  without  delay  present,  and  are 
the  aliasing  and  the  truncation  errors;  the  third  term  includes  the  error  from  the  cyclic 
delays  inherent  in  the  beamforming  process. 


APPENDIX  C 


ILLUSTRATION  OF  THE  METHOD 


A  computer  program  was  written  to  illustrate  the  theory.  A  set  of  random  numbers 
was  generated  using  a  CDC  Library  Function  Subroutine.*  This  set  was  then  digital  band¬ 
pass  filtered  for  C  =  W0/W  =  12  and  for  a  sampling  rate  SF  =  4(W0+W),  i.e.,  twice  the  usual 
Nyquist  rate.  The  resulting  set  of  numbers  is  taken  to  represent  a  redundant  sampling  of  a 
continuous  function  g(t)  bandlimited  to  the  interval  [Wj\,Wq+W].  The  set  can  then  be 
separated  into  two  related  interleaving  sets^V(t)  and  g^Ht),  either  set  completely  deter¬ 
mining  g(t)  (see  Fig.  Cl).  By  choosing  the  starting  point  for  g$(t),  the  set^2\f)  can  be 
“delayed”  from  o^V(f)  by  odd  multiples  of  1/3F.  For  the  case  shown  here  we  let  the  delay 
f0  between  the  two  sets  be  equal  to  13 fSF  =  A/2  =  (l/2W)/2. 

Both  sets  were  “bandpass  sampled”  at  a  rate  21 V  by  retaining  only  every  13th  value 
from  each  set,  thus  creating  two  sets  of  bandpass  sampled  data  delayed  from  each  other  by 
13 fSF  (see  Fig.  C2);  AT  was  equal  to  128,  and  Hanning  weighting  was  used  to  minimize  g(t) 
near  f=0  and  i-T.  Figure  C3  shows  a  plot  of  the  difference  between  the  original  Hanned  set 
values  before  delaying.  Note  the  large  difference  values,  indicating  that  (f)  and  g^\t) 
were  not  at  all  coincident  before  being  delayed. 


Fig.  Cl— A  set  of  randomly  generated  numbers  is  digitally  bandpass 
filtered  and  separated  into  two  related  interleaving  sets  f^(t)  *nd 
•*  shown  here.  Either  set  completely  determines  the  bandpass 
function  g{t). 


•Spectrum  analysis  of  the  sequence  before  filtering  resulted  in  a  Oat  power  spectrum,  indicating  that  the 
series  was  random  enough  for  our  purposes. 


24 


NRL  REPORT  73B9 


25 


Fig.  C2— The  set  g^\t)  shown  in  Fig.  Cl  is  delayed  from  by  an 

amount  tp  =  13 JSF.  Both  sets  are  then  “bandpass  sampled”  at  a  rate 
2V/  by  retaining  only  every  13th  value  from  each  set.  The  resultant  sets 
are  then  as  shown  here. 


N  =  120 


Fig.  C3— Amplitude  differences  between  the  original  Hanned  set  values 
g^y(t)  and  g^\t)  before  (top  curve)  and  after  (bottom  curve)  com¬ 
puting  the  advance  function  g^\t+to).  The  differences  are  large  before 
delaying  (top  curve).  After  delaying,  the  two  set  values  agree  closely 
(bottom  curve). 
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Equation  (16)  was  programmed,  for  f0  =  -13/SF,  and  run  on  the  g^Ht)  data  series. 
An  IDFT  of  the  resulting  Gji(fk')  was  computed,  and  the  values  of  the  “advanced”  func¬ 
tion  g^2)(f+|£0  I)  were  compared  to  those  of  A  plot  of  the  difference  function 

^V(t)-g^2)(t-Hfol)  is  shown  in  Fig.  C3.  It  is  seen  that  the  two  series  agree  closely,  having 
been  properly  delayed.  The  mean-square  value  of  the  “after-delay”  difference  function  is 
equal  to 

N 


el 


_1 

N 


n=nl=l 


Using  Parser'll ’s  theorem  which  equates  the  mean-square  value  of  a  function  to  the  sum  over 
all  frequencies  of  the  square  of  the  absolute  value  of  its  Fourier  spectrum,  and  with  the 
knowledge  that  the  spectrum  of  the  noise  is  approximately  flat  across  the  band,  we  com¬ 
pute  the  mean-square  error  per  spectral  estimate  to  be 


2c1 /N  =  5X10“6  =  —45  dB 


where  the  decibel  numbers  are  referenced  to  the  mean-square  value  of  the  original  time 
series  gj^t).  Using  Eq.  B16  we  find  that  for  our  choice  of  filter,  Hanning  smoothing 
window,  and  time  delay,  the  errors  from  aliasing,  truncation,  and  cycling  are  predicted  to 
be  less  than,  respectively,  -38  dB,  -65  dB.  and  -106  dB.  This  indicates,  then,  that  in  this 
case  almost  all  of  the  error  is  contributed  by  spectrum  aliasing.  There  is  almost  no  contribu¬ 
tion  to  the  error  from  the  delay  f0  -  A/2  (very  small)  or  the  truncation  (the  Hanning  window 
is  very  small  at  the  end  points  of  the  interval)  problem. 

The  same  two  original  data  sets  were  again  “bandpass  sampled”  and  properly  delayed, 
but  this  time  the  delay  was  very  large.  We  let  t0  =  (T/2)  -  (5A/2).  This  value  of  f0  is  close 
to  half  the  analysis  interval  T/2,  at  which  point  the  limiting  assumptions  for  negligibility  of 
g(t)  begin  to  break  down.  The  computed  variance  per  spectral  estimate  for  this  case  was  on 
the  order  of 


2e1/N  =  10-3  =  -22  dB. 

Again  using  Eq.  B16,  we  predict  the  errors  from  aliasing,  truncation,  and  cycling  to  be,  re¬ 
spectively,  —38  dB,  —65  dB,  and  —20  dB.  Clearly,  the  dominant  error  is  from  the  cyclic 
delay  term.  This  computed  error  agrees  closely  with  the  calculated  experimental  error  and 
shows  the  adequacy  of  the  analysis  of  App.  B  to  handle  this  extreme  delay  case. 
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